IMPROVING  LOW  ORDER,  LINEAR,  POSITIVE  SPATIAL 
QUADRATURES  FOR  THE  PARTIAL  CURRENT  NEUTRON  TRANSPORT 

METHOD 

THESIS 

John  M  Snyder,  Major,  USA 


AFIT/GNE/ENP/10M-08 


DEPARTMENT  OF  THE  AIR  FORCE 
AIR  UNIVERSITY 

Air  Force  Institute  of  Technology 

Wright-Patterson  Air  Force  Base,  Ohio 

APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


IV 


The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the 
official  policy  or  position  of  the  United  States  Air  Force,  Department  of  Defense, 
or  the  U.S.  Government. 


v 


AFIT/GNE /ENP/10M-M08 


IMPROVING  LOW  ORDER,  LINEAR,  POSITIVE  SPATIAL  QUADRATURES 
FOR  THE  PARTIAL  CURRENT  NEUTRON  TRANSPORT  METHOD 

THESIS 

Presented  to  the  Faculty 
Department  of  Engineering  Physics 
Graduate  School  of  Engineering  and  Management 
Air  Force  Institute  of  Technology 
Air  University 

Air  Education  and  Training  Command 
In  Partial  Fulfillment  of  the  Requirements  for  the 
Degree  of  Master  of  Science 

John  M  Snyder,  B.S. 

Major,  USA 

March  2010 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 


vi 


AFIT/GNE /ENP /10M-M08 


IMPROVING  LOW  ORDER,  LINEAR,  POSITIVE  SPATIAL  QUADRATURES 
FOR  THE  PARTIAL  CURRENT  NEUTRON  TRANSPORT  METHOD 


John  M  Snyder,  B.S. 
Major,  USA 


L0  to 

Kirk  A.  Mathews,  PhD  (Research  Advisor)  Date 

Professor  of  Nuclear  Engineering 

/Kt^c/ 2oLo 

f^&T-  LtCol  Kyle  A.  Novak,  PhD  (Member)  Date 

Asst.  Professor  of  Mathematics 


^L— ^  ~  10  3°/Q 

Cpt  Benjamin  R.  Kowash,  PhD  (Member)  Date 

Asst.  Professor  of  Nuclear  Engineering 


vii 


AFIT/GNE /ENP /10M-M08 


Abstract 


AFIT  researchers  have  developed  a  new  approach  to  solving  Discrete 
Ordinates  equations,  which  approximate  the  linear  Boltzmann  Transport 
Equation  (BTE).  The  usual  approach  is  von  Neumann  iteration  on  the  scattering 
source,  which  requires  repeated  sweeps  through  the  spatial-angular  grid. 
Acceptable  convergence  requires  complicated  and  expensive  acceleration  schemes. 
The  new  approach,  Partial-Current  Transport  (PCT)  with  Adaptive  Distribution 
Iteration,  eliminates  scattering  source  iteration  through  matrix  inversions  and  a 
reduced-size  global  linear  algebra  problem.  It  creates  the  needed  matrices  directly 
from  the  standard  spatial  quadratures  used  in  the  sweeping. 

Positivity,  linearity,  and  (higher-than- first-order)  accuracy  are  the  key 
desirable  qualities  with  all  Discrete  Ordinates  methods,  but  all  three,  according 
to  Lathrop  [8],  cannot  be  achieved  simultaneously.  If  a  high  order  accurate,  linear 
method  is  used,  it  can  produce  negative  fluxes.  Non-linear  methods  have  been 
developed  that  are  high-order  accurate  and  positive,  but  these  methods  are  not 
widely  accepted  because  the  BTE  is  itself  a  linear  equation.  Positive,  linear 
methods  are  available,  but  are  only  first-order  accurate.  The  latter  can  achieve 
needed  accuracy  by  using  optically-thin  cells,  but  with  Source  Iteration  (SI),  this 
requires  a  fine  grid  of  many  cells,  hence  large  computational  expense. 

My  new  approach  is  to  partition  an  optically  thick  cell  into  2N  identical 
sub-cells.  Each  sub-cell  is  optically  thin  enough  that  first-order  accurate  spatial 
quadrature  methods  are  sufficiently  accurate  as  well  as  being  linear  and  positive. 
The  needed  matrices  are  computed  as  before  for  a  (thinnest)  sub-cell.  My 
algorithm  combines  the  matrices  for  a  pair  of  sub-cells  to  get  the  matrices  for  a 
single  (merged)  sub-cell  twice  as  thick.  Merging  N  times  yields  the  matrices  for 
the  original  cell.  This  allows  PCT  to  solve  the  discrete  ordinates  equations  with 
linearity,  positivity,  and  sufficient  accuracy  without  the  high  computational  cost 
of  increasing  the  number  of  cells  by  a  factor  of  2N. 
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IMPROVING  LOW  ORDER,  LINEAR,  POSITIVE  SPATIAL  QUADRATURES 
FOR  THE  PARTIAL  CURRENT  NEUTRON  TRANSPORT  METHOD 

I.  INTRODUCTION 

Neutron  transport  has  been  a  building  block  of  nuclear  reaction  science 
since  the  discovery  of  fission  in  the  1930s  and  1940s,  both  for  calculations  of 
energy  output  and  for  personnel  safety.  The  neutron  is  the  necessary  catalyst  for 
nuclear  fission  and  the  output  of  import  of  both  nuclear  fission  and  fusion 
thereby  making  the  information  of  where  and  how  many  neutrons  there  are,  or 
neutron  flux,  a  primary  interest  to  nuclear  engineers. 

Solving  the  BTE  determines  the  amount  of  neutron  flux.  However,  the 
BTE  is  an  integro-differential  equation  that  is  not  directly  solvable  except  in  the 
simplest  cases;  therefore  it  requires  a  numerical  approximation.  Discrete 
Ordinates  (deterministic)  and  Monte  Carlo  (probabilistic)  methods  are  currently 
the  two  major  methods  for  approximating  the  BTE  for  neutron  transport. 

Discrete  ordinates  methods  have  been  popular  for  approximating  the  BTE  since 
the  early  days  of  nuclear  science  because  of  their  relative  computational  ease  [1]. 
For  efficiency  and  simplicity  linear  methods,  such  as  Diamond  Difference  (DD) 

[2],  which  is  2nd  order  accurate,  and  Linear  Discontinuous  (LD)  [6],  which  is  3rd 
order  accurate,  are  widely  used.  However,  unphysical  negative  fluxes  are 
artifacts  of  these  methods  in  many  common  scenarios.  Due  to  these  unphysical 
artifacts  discrete  ordinates  methods  are  rejected  by  many  (who  adopt  Monte 
Carlo  methods). 

Three  competing  issues  with  discrete  ordinates  methods  are  linearity, 
positivity,  and  accuracy  [8].  According  to  Lathrop,  positivity  can  only  be 
achieved  at  the  cost  of  accuracy  or  non-linearity.  Depending  on  the  optical 
thickness  and  the  scattering  and  absorption  properties  of  the  material  of  the 
problem,  discrete  ordinates  methods  can  produce  negative  fluxes  that  are 
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physically  nonsensical  and  computationally  time  consuming  to  correct.  Spatial 
distortions,  called  ray  effects,  also  arise  due  to  the  angular  quadrature  [7],  [13]. 
This  leads  to  the  use  of  more  complex  methods  in  order  to  minimize  the  amount 
of  negativity  [9]  or  to  confine  it  to  regions  of  low  concern.  In  order  to  increase 
the  accuracy  using  a  low  order,  linear  method,  the  spatial  mesh  must  be  refined 
to  a  point  that  becomes  computationally  expensive.  Figure  1  is  a  simple 
representation  of  an  angular  quadrature  producing  beams  or  “rays”  of  neutron 
flux. 


In  this  figure  pn  is  the  cosine  of  the  angle  between  the  direction  of  the  neutron 
path  and  the  x-axis  in  slab  geometry.  The  angular  quadrature  does  not  produce 
actual  “rays”  but  represents  the  number  of  points  used  to  approximate  values  on 
the  function,  like  the  number  of  points  that  are  used  to  approximate  an  integral 
using  the  Trapezoid  Method.  The  Sn  discretization  [1]  is  the  angular  quadrature 
approximation  that  I  use  throughout  this  project,  where  n  is  the  number  of 
directions  that  are  used  to  approximate  the  BTE  for  neutron  flux.  Another 
concern  with  this  type  of  angular  discretization  of  the  angular  flow  of  neutrons  is 
that  when  neutrons  collide  with  nuclei  there  is  a  continuous  distribution  of 
directions  into  which  they  can  scatter.  By  only  looking  at  certain  direction  of 
neutron  flow  this  limits  the  accuracy  of  the  approximation.  One  way  to  minimize 
the  difference  between  the  approximation  and  the  truth  is  to  increase  number  of 
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rays.  As  the  number  of  directions  approaches  infinity  the  approximation 
approaches  the  true  answer.  This  is  not  computationally  feasible.  However, 
increasing  the  number  of  ordinates  will,  in  theory,  bring  the  approximation  closer 
to  the  truth. 

In  recent  years  many  new  methods,  as  well  as  enhancements  to  old  methods, 
have  been  investigated  to  minimize  the  negativity  and  increase  the  accuracy  of 
the  answer.  Mathews  et  al  [15]  have  developed  adaptive  split-cell  characteristic 
methods  of  3rd  or  4th  order  that  preserve  positivity.  However,  these  methods, 
which  are  non-linear,  increase  the  computational  costs,  particularly  with  the  new 
partial-current  transport  method  for  Sn  [5]. 

The  Adaptive  Partial-Current  Discrete  Ordinates  Radiation  Transport  with 
Distribution  Iteration  (DI)  [5],  [14],  [17],  [19]  solves  many  of  these  problems. 
However,  taking  advantage  of  the  linearity  of  the  BTE,  attaining  the  accuracy  of 
the  higher  order  methods,  and  not  introducing  the  negative  fluxes  associated  with 
many  spatial  quadratures  requires  a  method  for  DI  that  refines  the  spatial  mesh 
but  does  not  increase  the  computational  cost.  In  this  research  I  develop  an 
approach  to  solve  the  within-cell  transport  at  negligible  extra  computational  cost 
for  DI  using  the  Sn  angular  quadratures  using  low  order,  linear,  unconditionally 
positive  spatial  quadrature  methods  in  a  way  that  achieves  the  accuracy  of  higher 
order  methods  without  the  negative  flux  artifacts  of  those  high  order  methods.  I 

do  this  by  using  creating  a  very  optically  thin  sub-cell,  —the  thickness  of  the 

2n 

original  cell,  where  the  low  order  quadratures  are  very  accurate,  and  then  use  a 
merging  scheme  to  combine  the  cells  back  to  the  original  cell  thickness  while 
retaining  the  accuracy  of  the  thin  sub-cell. 

I. A.  Background 

The  issues  that  have  plagued  the  discrete  ordinates  methods,  negative  fluxes 
and  accuracy,  have  given  rise  to  many  attempts  to  minimize  the  presence  of  these 
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artifacts  and  to  improve  the  accuracy  of  the  solution.  The  preferred  method  of 
solving  the  neutron  distribution  in  a  material  is  through  the  SI  method.  SI  is  a 
Fixed  Point  Iteration  on  the  neutron  source.  This  method  uses  a  combination  of 
an  angular  quadrature  plus  a  spatial  quadrature,  uses  a  guess  for  the  source,  and 
then  walks  from  one  side  of  the  problem,  computing  the  neutron  flux  in  each  cell, 
to  the  other  side  of  the  problem.  With  the  neutron  flux,  the  balance  equation 
can  then  be  solved  for  the  source.  The  newly  computed  source  is  then  compared 
to  the  original  guess.  This  is  repeated  until  the  source  is  within  a  given  tolerance 
of  the  last  iteration.  As  the  problem  becomes  optically  this  or  the  material 
approaches  a  pure  scattering  source,  the  number  of  iterations  required  to  achieve 
a  convergence  tolerance  increases.  Thus,  attempts  to  improve  discrete  ordinates 
have  looked  at  ways  to  enhance  the  SI  method. 

If  the  spatial  quadratures  are  limited  to  low  order,  linear,  and  conditionally 
positive  the  spatial  mesh  must  be  refined.  This  requires  SI  to  iterate  through 
more  cells,  slowing  the  convergence.  A  method  that  maintains  all  three  of  the 
desirable  qualities  has  not  been  developed  where  the  benefits  (increase  in 
accuracy)  don’t  outweigh  the  cost  (increase  the  number  of  calculations  slowing 
the  convergence).  Different  basis  functions  for  the  spatial  and  angular 
quadrature  have  been  developed.  In  this  project  I  do  not  look  in  depth  at  SI,  but 
instead  my  method  improves  the  performance  of  PCT  with  Adaptive  DI  by 
improving  the  fidelity  of  the  spatial  model  (for  a  given  Sn  model)  of  transport 
within  each  cell. 

PCT  with  Adaptive  DI  looks  at  the  problem  in  a  different  light.  Instead  of 
iteratively  sweeping  through  the  mesh  until  the  solution  converges  to  an  analytic 
solution,  DI  solves  the  within-cell  transport  by  creating  a  transport  matrix  for 
the  entire  cell  that  is  based  on  the  characteristics  of  the  medium.  This  matrix 
describes  all  the  transport,  including  scatter  and  emission  through  the  cell.  With 
the  internal  transport  of  each  cell  correct,  the  inter-cell  transport  can  then  be 
solved  using  the  partial  currents  at  each  face.  This  still  requires  some  iteration, 
but  far  fewer  in  most  cases  than  comparable  SI  schemes.  DI  shows  promise  but 
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still  may  fail  to  perform  effectively  in  some  problems,  specifically  2-D  and  3-D, 
when  spatial  quadratures  that-  are  not  unconditionally  positive  are  used  [12]. 

This  leads  to  the  necessity  of  providing  DI  with  a  way  to  increase  the  accuracy  of 
positive,  linear  spatial  quadratures  to  the  same  level  as  higher  order  methods 
that  produce  unphysical  results. 


I.A.l  Linearity 

The  BTE  for  neutrons  is  a  linear  equation.  Therefore  any  numerical 
approximation  to  the  linear  equation  should  itself  be  linear.  There  are  non-linear 
approximations  to  the  BTE  that  have  been  shown  to  be  3rd  or  4th  order  accurate. 
However,  these  have  not  been  widely  accepted  because  in  certain  geometries  they 
do  not  reach  the  appropriate  diffusion  limit  [11],  Also,  because  these  quadratures 
are  non-linear,  after  each  iteration,  some  parameters  of  the  method  change 
because  the  flux  changes,  so  all  the  spatial  quadrature  coefficients  must  be 
recalculated  and  new  matrices  (for  DI)  are  needed  and  are  different  in  every  cell. 
This  is  disastrous  in  both  storage  (RAM)  and  time.  Therefore,  methods  that  are 
linear  (i.e.  DD,  LD,  Linear  Nodal  (LN),  and  Linear  Characteristic  (LC))  are  more 
accepted  and  generally  more  desirable. 


I.A.2  Positivity 

As  pointed  out  by  Lathrop  [8]  the  primary  trade  off  for  discrete  ordinates 
methods  comes  in  part  due  to  the  physical  requirement  of  the  solution  (i.e.  the 
neutron  flux,  or  dose)  to  be  positive.  Negative  fluxes  arise  because  of  a)  Too 
thick  cells  when  using  conditionally  positive  spatial  quadratures  (e.g.  DD  in  slab 
geometry),  b)  Discontinuous  boundary  conditions  with  2-D  or  3-D  and  methods 
with  negative  coefficients  (e.g.  DD  on  rectangles),  and  c)  Truncated  Spherical 
Harmonics  /  Legendre  expansions  of  anisotropic  scatter  combined  with  an 
unfortunate  choice  of  angular  quadrature  set.  Many  scientists  and  engineers  will 
not  blithely  accept  the  condition  that  there  is  a  negative  flux  of  neutrons  on  the 
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basis  that  the  negative  flux  in  question  is  present  in  only  a  portion  of  the  solution 
space.  This  casts  serious  doubts  onto  the  validity  of  the  method  as  a  whole.  The 
fix-ups  that  have  been  developed  look  for  negative  fluxes  and  then  adapt  the 
method  to  that  space,  forcing  the  solution  to  zero  or  a  positive  value,  then 
rebalance  the  equation.  By  rebalancing  the  equation,  the  fix  up  ensures  particle 
conservation,  but  it  also  reduces  the  method  to  1st  order  accuracy  wherever  the 
fix  up  is  applied. 


I.  A. 3  Accuracy 

The  final  issue  of  importance  for  discrete  ordinates  methods  is  that  of 
accuracy.  Again  pointing  to  Lathrop’s  discussion  of  the  tradeoff  between 
positivity  and  accuracy,  the  2nd  order  and  higher  methods  that  are  linear  are  at 
best  only  conditionally  positive.  DD  is  not  even  conditionally  positive  in  2-D  and 
3-D.  The  fix-ups  that  have  been  employed  to  overcome  the  negativity  in  some 
quadrature  methods  are  generally  only  1st  order  accurate.  By  using  these  fix-ups 
the  accuracy  gained  by  the  higher  order  conditionally  positive  methods  are 
negated  by  the  low  order  of  the  fix-up.  This  shows  that  of  the  three  desirable 
characteristics  of  discrete  ordinates  methods,  it  is  possible  to  have  only  two  out 
of  three.  You  can  have  high  order  accuracy  and  positivity,  but  not  linearity.  Or 
you  can  choose  high  order  accuracy  and  linearity,  but  not  guarantee  positivity. 
And  finally,  you  can  have  linearity  and  positivity,  but  attain  only  1st  order 
accuracy. 


I.A.4  Distribution  Iteration 

DI  was  developed  by  Mathews,  Prins,  Wager,  and  Dishaw  [5], [14],  [17],  [19], 
with  the  idea  that  the  within-cell  transport  could  be  represented  by  a  set  of 
coefficient  matrices  that  act  on  the  incoming  flux  and  contain  within  them  all  the 
information  regarding  the  transmission,  reflection,  scattering  and  escape  of 
neutrons  through,  within,  and  out  of  the  medium.  The  matrices  are  set  up  using 
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a  set  of  balance  equations  dependant  on  the  properties  of  the  medium,  thickness 
of  the  cell,  and  the  spatial  quadrature  that  is  being  used.  The  matrices  can  be 
set  up  once  and  then  they  do  not  change  for  that  medium.  In  this  project  I  will 
use  the  Sn  angular  quadrature  set  to  approximate  the  BTE  within  each  cell.  For 
example,  Figure  1  shows  an  S4  representation,  in  which  there  are  four  discrete 
angles,  pi,  p2,  p3,  and  p4. 


I.B.  Motivation 

A  method  needs  to  be  developed  that  enables  a  linear,  positive  spatial 
quadrature,  that  is  low  order  accurate  to  attain  the  accuracy  of  the  higher  order 
quadratures  that  violate  linearity  or  positivity.  This  will  bypass  the  possible 
failures  that  may  arise  in  DI  from  quadratures  that  are  non-linear  or  that 
produce  negative  neutron  flux. 


I.C.  Objectives 

The  objectives  of  this  research  are  to:  1)  increase  the  accuracy  of  positive, 
linear  spatial  quadratures  to  that  of  higher  order  methods,  2)  maintain  numerical 
stability  in  the  matrices  that  it  produces,  and  3)  produce  a  method  that  can  be 
integrated  into  existing  DI  codes. 

1.  Increased  Accuracy:  Based  on  the  desire  to  have  a  linear  and  positive 
spatial  quadrature,  and  due  to  the  fact  that  such  methods  are  only  1st  order 
accurate,  this  is  the  key  objective  of  this  project.  The  results  will  be  compared  to 
a  benchmark  to  validate  that  the  method  can  produce  similar  accuracy  of  higher 
order  methods.  This  method  will  be  compared  to  a  benchmark  (SI)  code  that- 
has  been  developed  and  improved  over  the  last  15  years  [18].  This  will  insure 
that  this  method  can  achieve  the  accuracy  of  higher  order  methods. 

2.  Numerical  Stability:  In  order  for  this  method  to  reach  the  accuracy  of 
higher  order  methods,  matrix  inversions  are  required  to  emulate  the  mesh 
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refinement  required  by  SI  methods.  This  will  require  that  the  matrices  be  well- 
conditioned. 

3.  Integration  into  DI  Code :  Upon  successful  testing,  my  method  will 
produce  the  matrices  that  are  used  in  slab  geometry  by  DI  codes.  This  method 
will  improve  the  performance  of  1st  order,  linear,  positive  spatial  quadratures  in 
the  existing  DI  codes. 


I.D.  Statement  of  the  Problem 

Develop  a  method  that  can  be  incorporated  into  existing  slab  geometry  DI 
codes  that  implements  linear,  unconditionally  positive,  hence  only  1st  order 
accurate,  spatial  quadratures  for  the  discrete  ordinates  approximation  to  the 
BTE  that  achieves  the  accuracy  of  high  order  methods  that  are  not  positive  or 
not  linear. 


I.E.  Scope  and  Limitations 

Method: 

1-D  Cartesian  (Slab)  geometry 
1  Energy  group 
Testing: 

Subcritical  steady  state  (not  eigenvalue  problems) 

C  <  1 

Demonstrate  Step  and  SC,  benchmarking  each  against  SI 
Compare  to  DD  and  EC 


I.F.  Approach 

The  first  step  is  to  set  up  the  equations  for  the  partial  currents  and  fluxes  for 
the  original-sized  cell  based  on  the  spatial  and  angular  quadratures.  Once  the 
cells  are  split  the  next  step  is  to  solve  the  simultaneous  equations  for  the  two 
cells  in  order  to  eliminate  the  currents  and  fluxes  between  the  two  adjacent  cells, 
thus  combining  the  two  smaller  cells  into  a  cell  twice  as  big.  This  gives  the 


equations  for  the  combined  cells,  and  is  repeated  until  the  sub-cells  are  combined 
to  the  original-sized  cell. 

The  next  step  is  to  use  the  Step  and  SC  spatial  quadratures  (linear,  positive, 
and  1st  order)  and  test  and  verify  the  currents  and  fluxes  produced  by  this 
method  against  the  same  quadratures  using  an  SI  code. 

Finally  the  method  is  compared  to  a  benchmark  (Exponential  Characteristic 
(EC)  in  the  SI  code)  in  several  different  problem  scenarios  to  show  that  the  low 
order,  positive,  linear  spatial  quadratures  can  provide  the  accuracy  of  the  higher 
order  method.  The  problem  scenarios  use  characteristics  of  cell  width  (based  on 
the  mean  free  path  of  the  neutrons  in  the  given  material),  scattering  ratios, 
internal  neutron  source,  and  external  neutron  illumination  to  compute  the  partial 
currents  out  each  face  of  the  original  cell  and  the  average  scalar  flux  in  the  cell. 
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II.  THEORY 


As  all  discrete  ordinates  methods,  and  neutron  transport  theory  in  general, 
the  object  of  the  approximation  is  to  find  the  flux  of  neutrons  at  a  particular 
location  or  throughout  the  space  of  interest.  The  BTE  is  the  basis  for  all  neutron 
transport  and  in  its  full  form  describes  the  motion  of  neutrons  based  on  their 
energy,  location,  and  the  scattering  and  absorption  properties  of  the  material 
through  which  they  are  moving.  Without  making  any  simplifying  assumptions 
regarding  energies  or  discretizing  space  and  angle  we  assume  the  BTE  for 
neutron  transport  to  be: 

+  ^  '  V  +  E,t)  =  q(r,Cl,E,t).  (1) 

V  ot 

For  this  work  the  BTE  has  been  reduced,  in  discrete  ordinates,  to  a  time- 
independent,  mono-energetic  cell  balance  equation  in  the  following  form: 

jOllt  +  G  ^ x  yy  ^  ^  _|_  SeXtAX  +  jm,  (2) 


where  cr  is  the  total  cross  section  of  the  material,  oy  is  the  scattering  cross 
section  and  Ax  is  the  thickness  of  the  cell.  Figure  2  is  a  representation  of  a  cell 
with  the  transport  of  neutrons  into,  through  and  out  of  the  cell  in  all  directions 
(//) ,  and  where: 


A  ¥t(m) 

"A  VL{n) 


//  >  0, 

A  <  0, 


(3) 


yL(lu)\s  the  flux  at  the  left  side  of  the  cell;  y/R{/u)is  the  flux  at  the  left  side  of 
the  cell;  y/ 4(a)  is  the  average  flux  in  the  cell  (spatially  averaged);  Sext(ju )  is  the 

emission  source  inside  the  cell,  assumed  to  be  spatially  uniform  and  isotropic;  and 

l 

(p ^  =  J  y/A(ju)djU  is  the  average  scalar  flux  in  the  cell. 

-l 
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Figure  2:  Within-Cell  Transport 


Neutron  transport  does  not  include  anti-neutrons  or  neutron  holes  (unlike 
electron  transport  in  transistors),  so  the  actual  flux  is  non- negative: 

V(x,n)  >  0. 


(4) 


A  numerical  method  that  ensures  that  its  solutions  have  this  property  is  said  to 
be  a  positive  method.  In  discrete  ordinates  calculations,  most  spatial  quadrature 
methods  are  not  positive. 


II. A.  Angular  Quadratures 

There  are  many  ways  to  discretize  the  angular  space  of  neutron  transport 
problems.  For  this  project  I  use  the  Sn,  or  discrete  ordinates,  angular  quadrature 
approximations,  specifically  the  single-range  Gauss-Legendre  angular  quadrature 
sets  in  slab.  See,  for  example,  Lewis  and  Miller  [7].  The  approximation 
approaches  the  analytical  solution  to  the  transport  equation  as  the  number  of 
ordinates  approaches  infinity.  This  is  not  computationally  possible.  Although  I 
explored  the  performance  of  my  methods  using  S2  through  S8,  the  testing 
reported  here  used  S8,  being  the  most  challenging  of  these. 
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II. B.  Spatial  Quadratures 

As  the  problem  statement  declares,  the  crux  of  the  problem  is  to  develop  a 
method  that  allows  a  low  order  spatial  quadrature  that  is  linear  and 
unconditionally  positive  to  achieve  the  accuracy  of  a  high-order  quadrature  that 
may  be  either  linear  or  positive,  but  not  both.  While  many  spatial  quadratures 
have  been  proposed,  I  used  four  spatial  quadratures  in  this  work:  DD,  Step,  Step 
Characteristic  (SC),  and  EC.  Table  1  summarizes  the  four  spatial  quadratures 
on  which  I  focus. 


Table  1:  Spatial  Quadrature  Comparison 


Name 

Abbreviation 

Order* 

Positivity 

Linearity 

Description 

Step 

ST 

1 

X 

X 

Constant  discontinuous 
flux  approx. 

Step 

Characteristic 

SC 

1 

X 

X 

Constant  discontinuous 
scattering  source  approx 

Diamond 

Difference 

DD 

2 

X 

Divided  difference  approx. 

Exponential 

Characteristic 

EC 

4 

X 

*  error  proportional  to 

((TAX)71 


Of  the  low-order  methods,  DD  is  the  most  accurate,  having  second  order 
accuracy,  but  is  only  conditionally  positive.  The  condition  that  gives  DD  the 
negativity  is  the  thickness  of  the  cell  in  which  it  is  used.  This  condition  takes  the 
DD  method  out  of  consideration  for  this  method.  Interestingly,  even  though  DD 
has  this  conditional  positivity,  it  is  still  the  preferred  choice  of  some  members  of 
the  transport  community. 

Exponential  Characteristic  is  a  method  that  was  developed  for  slab  geometry 
at  the  Air  Force  Institute  of  Technology  [15].  This  method  uses  a  characteristic 
integration  of  the  Boltzmann  transport  equation  with  an  exponential  function  as 
the  assumed  form  of  the  source  distribution,  continuous  across  each  spatial  cell. 
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This  leaves  two  methods,  both  of  which  are  only  first  order  accurate,  but  are 
both  unconditionally  positive,  Step  and  Step  Characteristic  (SC).  Step  is  a  linear 
discontinuous  approximation.  The  Step  approximation  is  similar  to  a  backward 
Euler  approximation.  Step  makes  the  auxiliary  assumption  that  the  average  flux 
in  the  cell  is  constant  and  continuous  at  the  outflow  of  the  cell  but  discontinuous 
at  the  inflow  of  the  cell.  SC  makes  a  similar  assumption  of  a  constant  scattering 
source  that  is  that  is  continuous  at  the  outflow  of  the  cell  but  discontinuous  at 
the  inflow  of  the  cell.  Of  these  two  methods  step  is  the  most  straightforward, 
but  neither  one  is  complicated  or  expensive  to  calculate.  Based  on  this  summary 
of  quadratures  and  being  consistent  with  the  objectives  of  this  project,  it  is  clear 
that  the  only  spatial  quadratures,  of  the  four  examined  here,  that  meet  the 
criteria  are  Step  and  SC. 


II.  C.  Within- Cell  Transport 

The  within-cell  transport  problem  is  at  the  heart  of  the  PCT  method:  It 
solves  the  within-cell  transport  and  collapses  that  into  a  global  partial  current 
that  is  used  to  determine  the  partial  currents  out  of  each  face  of  a  cell  based 
upon  the  partial  currents  into  each  face  of  the  cell  and  the  internal  source. 


II.C.l  Calculating  the  Within- Cell  Transport  Matrices 
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While  much  of  the  following  equations  were  developed  for  DI  by  Mathews  and 
Dishaw  [5],  [14],  they  are  critical  to  this  project.  Therefore,  an  explanation  of  the 
terms  and  algebraic  manipulation  of  the  matrices  involved  is  advantageous.  The 
solution  to  the  method  can  be  observed  as  either  the  currents  coming  out  of  the 
cell  or  as  the  average  flux  in  the  cell.  Because  the  flux  in  the  cell  generated  a 
scattering  source  which  in  turn  creates  a  flux  or  current  the  two  systems  are 
closely  related.  I  will  develop  both  of  the  system  of  equations  as  they  relate  to 
this  method.  The  basic  algorithm  for  determining  the  currents  out  of  each  of  the 
cell  faces  and  the  average  flux  in  the  cell  is  as  follows: 


Starting  with  one  cell  as  in  Figure  2,  where,  for  S4  (with  -1  <  <  ..  <  <  1 


J°Ut  = 


f  r  \ 

'  jLll 

■L 

jR 

JjU?j 

jR 

vV4/ 


(5) 


and 


jm  = 


f  jR  A 

iR 

■l/u2 

i/uZ 


(6) 


The  average  flux  in  the  cell  can  be  represented  as 

y/  =  +  K^’^S  ext'  +  K^’^S5, 


(7) 


where  vp  is  the  average  flux  in  the  cell,  Sext  is  the  source  due  to  external  emission 
within  the  cell,  and  Ss  is  the  source  due  to  scattering  within  the  cell,  all  of  which 
are  vectors  of  size  n  (Sn),  based  on  the  angular  quadrature  used.  Kv/,m 
represents  the  flux,  approximated  by  a  specific  spatial  quadrature,  that  is  due  to 
current  in.  Likewise,  KV/,E  and  represent  the  flux  due  to  the  external 
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source,  and  the  flux  from  the  scatter  within  the  cell,  respectively.  The  K 
matrices,  developed  in  Appendices  A  and  B  for  Step  and  SC,  respectively,  are 
square,  diagonal  matrices  (n  x  n). 

The  scatter  source  is  a  scattering  matrix  that  operates  on  the  flux 

Ss  =  Z^y  (8) 

Substituting  (8)  into  (7)  gives  the  flux  in  terms  of  itself, 

)//  =  K^infn  +  K^ESext  +  K^Z5’*>  (9) 

Solving  for  v|/  using  matrix  notation  gives 

Iv|/  -  K^Z5’^  =  K^’*nj*n  +  K¥’ESext  (10) 

V)/  =  (I  -  K^Z^r^K^j™  +  Kv'-ESext  )  (11) 

The  inverse  operator  in  this  equation  accounts  for  all  scatters  that  occur 
throughout  the  cell.  Let 

L  =  (I  -  Kv',s'Zs,v/)~l  (12) 

So  that 

V  =  LK^’injin  +  LK^S63*  (13) 

Now  the  flux  that  is  generated  by  the  external  source  and  the  flux  that  is 
generated  directly  from  current  coming  into  the  cell  can  be  separated.  This  gives 
m  matrices  for  the  flux 

mys,in  =  LKy,,in  (14) 

and 

m r,ext  =lkV,,e  '  (i5) 

Substituting  (14)  and  (15)  into  (13)  gives  the  solution  for  the  flux  in  the  cell  from 
current  coming  into  the  cell  and  external  emission  source: 

V  =  fn  +  m^extSext  (16) 
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Looking  at  the  contributions  to  currents  out  of  the  cell  that  come  from  all 
sources  is  similar  to  the  contributions  to  the  flux: 


_  j^rO,in j in  j ^0,E^ext  _j_  (1^) 

Substituting  (8)  in  to  (17)  gives 

j  °ut  =  K°’Jjin  +  }^0,Esext.  +  }^0,S^S,<y^  (18) 


Now  by  substituting  (13)  into  (18)  gives  the  total  current  out  of  the  cell  based 
current  into  the  cell,  current  from  the  external  source  and  the  current  from  the 
scatters: 


jO«£  _ 


-^rO,in  jin 


+  K°’E  g  ext 


(LK^’ 


in  jzn 


+ 


LK^S' 


j  ext\ 


(19) 


Expanding  and  rearranging  (19)  gives 


•  out 


_  (K^’*”  +  )  j*n  _|_  (K°,e  + 


kO^S'Es^LKv/’E)s 


ext 


(20) 


Collecting  the  terms  that  act  on  jmand  Sext  gives  the  m  matrices  that  act  on 
jmand  Sext  in  order  to  produce 

~»out  _  rnOut,in~»in  .  out,extQext  -i  \ 

J  —  HI  J  -h  III  O 


where  the  matrix  of  coefficients  that  acts  directly  on  jm  represent  the  m 


out, in  , 


m 


out, in  j ^0,in  _j_ 


(22) 


and  the  matrix  of  coefficients  that  act  on  the  external  emission  source  and  source 
produce  by  scattering  is: 

m°ut,ext  =  kO,E  +  kO,5sS,^lk^,E  (23) 


The  m  matrices  can  also  be  represented  as  the  transmittance  and  reflectance  of 
the  neutrons  that  pass  through  the  cell  in  the  case  of  jm  ,  or  the  escape  of 
neutrons  that  are  emitted  in  the  cell  in  the  case  of  Sext ,  producing  either  current 
out  of  the  cell  through  each  of  the  cell  faces,  to  include  any  scattering  that  may 
happen,  or  the  average  flux  in  the  cell.  The  m  matrices  for  the  currents  out  of 
the  cell  are 
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(24) 


and 


m 


out, in  _ 


r  T  r~ 

j 

R  + 

V  J 


R 

j 

rji  +  + 

j  J 


out, ext 


( 

E . 

j 

E  +~ 
V  3 


V+J 


(25) 


In  these  matrices  the  superscripts  indicate  the  final  direction  of  flow  followed 
by  the  initial  direction  of  flow  and  the  subscripts  specify  the  source  of  the  flow, 
either  from  current,  j  ,  or  flux,  y/  .  T  corresponds  to  a  transmission  of  neutrons 
through  the  cell,  R  to  a  reflection  in  a  cell,  and  E  to  the  emission  of  neutrons 
from  the  external  source,  Sext .  For  example,  T  1  is  the  transmission  of 
neutrons  from  flux  that  was  initially  going  in  the  positive  direction  but  after 
some  number  of  scatters  is  now  going  in  the  negative  direction. 

In  equation  (24)  IV  the  transmission  of  jm  that  enters  the  cell  through  the 
right  boundary,  moving  in  the  negative  x  direction  (one  of  the  negative  p 
ordinates),  is  scattered  zero  or  more  times,  and  exits  through  the  left  boundary  in 
the  negative  direction  (continuing  in  one  of  the  negative  p  ordinates).  R  r+  is 

the  reflection  of  f  "  that  enters  the  cell  at  the  left  boundary  moving  in  the 
positive  x  direction  (one  of  the  positive  p  ordinates),  is  scattered  one  to  many 
times  and  leaves  through  the  left  boundary  moving  in  the  negative  x  direction 
(one  of  the  negative  p  ordinates).  Similarly,  R  is  the  reflection  of  jm  that 

enters  in  at  the  right  boundary  moving  in  the  negative  direction  (one  of  the 
negative  p  ordinates),  scatter  one  to  many  times  and  leaves  the  cell  through  the 
right  boundary  moving  in  the  positive  x  direction  (one  of  the  positive  p 
ordinates).  Finally,  TV++  is  the  transmittance  of  jm  through  the  cell  with  zero 
or  more  scatters,  entering  the  cell  through  the  left  boundary  moving  in  the 
positive  x  direction  (one  of  the  positive  p  ordinates)  and  leaving  the  cell  through 
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the  right  boundary  continuing  in  the  positive  x  direction  (continuing  in  one  of 
the  positive  p  ordinates).  I  have  used  the  mou  ,in  matrix  as  an  example. 

However,  the  other  matrices,  mout’ext )  m^’m  .  and  ;  are  similar. 

Special  care  must  be  taken  to  ensure  that  the  correct  sub-matrices,  T,  R, 
and  E,  are  used  in  these  equations  because  they  are  unique.  This  complicates  the 
algebra  but  once  the  system  is  set  up  it  can  be  repeated  easily.  Substituting  the 
sub-matrices  for  the  m  matrices,  we  can  now  represent  the  current  as  follows: 


T  ~~rR  +  R  ~+]+L  +  E  -S^’-  +  E  ~+Sext’+ 

3  J  3  J  3  3 

R  +~rR  +  T  ++]+L  +  E  +~sext~  +  E  ++Sext-+ 

3  J  3  J  3  3 


(26) 

(27) 


where  now  each  of  the  j  and  S'"'  vectors  is  of  length  n/2  and  sub- matrices  are  of 
size  (n/2  x  n/2). 

The  flux  can  also  be  broken  down  into  it  individual  sub-matrices  acting  on 
the  current  into  the  cell  and  the  external  emission  source: 


w~  =  T  “~r  +  R  “ +\+L  +  E  ~“Se 

T  i//  °  y/  °  y/ 


E  “ +Se 

v 


y/+  =  R  +  T  ++yL  +  E  1 

r  y/  °  y/  °  y/ 


-Qext,—  _j_  ^  ++gex^.+ 


(28) 

(29) 


Algorithm  1  is  the  basic  process  that  DI  uses  to  solve  the  partial  currents. 
This  is  done  on  a  fixed  size  cell  which  only  allows  the  spatial  quadrature  to 
attain  the  accuracy  for  that  optical  thickness.  My  method  of  merging  optically 
thin  cells,  demonstrated  in  the  next  section,  can  be  used  by  DI  to  solve  the  same 
problem  but  allows  the  first-order,  positive  spatial  quadratures  to  achieve  a  much 
higher  accuracy. 
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Algorithm  1:  Basic  Conceptual  Algorithm 


Initialize : 

1.  Use  problem  data  and  the  choice  of  basis  functions  to  generate  the  K 
matrices . 

2.  Solve  for  the  cell  m  matrix. 

3.  Use  emission  sources  and  within-cell  transport  to  evaluate  outflow 
currents  due  to  emission  sources 

4.  Initialize  the  inflow  current  distributions  (e.g.  uniform  and 
isotropic) 

5.  Solve  the  for  the  partial  currents  out  each  face  and  average  flux 
within  the  cell 


II.  C.2  Merging  Cells 

The  previous  section  characterizes  the  transport  of  neutrons,  in  an  angular 
discretization  of  the  cell,  based  on  an  arbitrary  spatial  quadrature,  through  one 
cell.  The  accuracy  of  the  computations  is  implicit  in  the  K  matrices  that  are 
dependent  upon  cell  thickness,  Ax  .  As  discussed  in  the  background  section,  in 
order  for  a  linear  and  positive  spatial  quadrature  to  attain  high  accuracy  the  cell 
must  be  sub-divided  until  it  is  optically  thin,  in  order  for  the  first-order,  positive, 
linear  methods  to  be  accurate.  In  the  SI  method,  this  is  done  by  refining  the 
mesh  and  sweeping  through  more  cells.  However,  as  pointed  out,  this  slows  the 
computation  time  and  when  extended  to  2-  and  3-D  geometries,  increases  the 
cost  of  the  calculation  faster  than  the  benefits.  As  will  be  shown  this  project  will 
produce  a  method  that  calculates  the  m  matrices  for  the  smallest  sub-divided  cell 
and  merges  them  matrices  back  to  the  original  size  cell. 

A  simple  illustration,  shown  in  Figures  4  through  11,  will  assist  in 
understanding  the  method.  Using  n  =  3  the  original  cell  is  divided  into  23  (8) 
sub-cells: 
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< -  Ax  - > 


Figure  4:  Original  Cell 


DO  I  =  n,  1,  -1 

i  =  3:  Merge  two  of  the  smallest  sub-cells 


Rearrange  terms  and  collect  T,  R,  and  E  sub-matrices  into  new  m 


matrices.  This  now  represents  a  set  of  sub-cells  twice  as  big. 


Figure  7:  Merge  1  Result 
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n  =  2:  Merge  two  of  the  smallest  sub-cells 


Figure  8:  Merge  #2 


Rearrange  terms  and  collect  T,  R,  and  E  sub-matrices  into  new  m 


matrices. 


Figure  9:  Merge  2  Result 


n  =  1:  Merge  two  of  the  smallest  sub-cells 


Figure  10:  Merge  #3 


n  =  0:  Rearrange  terms  and  collect  T,  R,  and  E  sub-matrices  into  new  m 
matrices. 


< -  Ax  - > 


Figure  11:  Merge  3  Result 
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Comparing  Figure  11  to  Figure  4,  after  3  merges  the  cell  is  back  to  its  original 
thickness.  However,  it  now  has  the  same  accuracy  for  the  solution  to  the  current 
out  and  the  flux  in  the  original  cell  as  if  it  was  1  / 8th  the  thickness.  On  a  scale  of 
23  this  does  not  appear  to  be  that  important.  However,  if  the  thickness  required 
for  accuracy  is  210,  SI  will  be  required  to  compute  the  current  and  flux  through 
1,024  sub-cells,  whereas  the  method  for  DI  will  compute  only  10  merges  and 
merges  are  only  done  once  for  one  cell  of  a  specific  size  and  composition 
(material).  It  will  be  shown  later  that  using  Step  in  SI  for  a  particular  problem 
required  219  (524,288)  sub-cells  in  order  to  attain  the  same  accuracy  as  using  EC 
in  SI.  For  this  method  that  requires  only  19  merges,  this  constitutes  a  savings  of 
computation  time  of  many  orders  of  magnitude. 

Using  my  method  I  am  able  to  accomplish  this  without  an  expensive  increase 
in  the  computation  time,  since  all  of  these  calculations  are  done  once  at  the 
beginning  of  the  problem  and  instead  of  having  to  sweep  through  2n  sub-cells,  as 
in  SI,  my  method  only  has  to  calculate  n  merges  in  order  to  achieve  the  same 
level  of  accuracy. 


II.  C.  2.1  The  Merging  Process 

Figure  12  shows  simple  example  of  dividing  the  original  cell  into  two  equal 
cells  of  half  the  original  size.  Even  though  the  picture  is  simple  the  algebraic 
manipulation  of  the  matrices  is  not.  In  this  section  I  will  explain  each  step  of  the 
merging  process  of  two  identical  sub-cells  that  can  be  done  repeatedly  until  the 
cell  is  back  to  the  original  thickness. 


rL , 

rc 

rR 

j-R 

(rc 

.  3  L 

Figure  12:  Currents  in  Split  Cells 
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Using  the  same  scheme  as  before,  the  m  matrices  can  be  formed  for  the 
smallest  sub-cell.  Since  this  cell  is  comprised  of  the  same  material  as  before  the 
computation  of  the  K  and  m  matrices  can  be  done  using  the  same  algorithm. 
Using  Figure  11  as  an  example,  where  n  =  1,  the  current  out  of  this  new  system 
of  cells  can  be  represented  the  same  as  in  (26)  and  (27). 

rL  =  T  ~rc  +  r  r+]+L  +  e  “ Sext~  +  E  ~+sext’+ 

J  3  J  3  J  3  3 

]+c  =  R  +~YC  +  T ++j+i  +  E  +~sext~  +  E  ++Sext+ 

J  3  J  3  J  3  3 

j ~c  =  T.  +  R  r+j+c  +  E  “Se:rf“  +  E  -+sext’+ 

J  3  J  3  J  3  3 

•+i?  H —  •  — R  rji  H — f  »+(7  -|-  J)  ^ ^  +  +  g63l£.+ 

Because  the  material  is  the  same  for  both  of  the  sub-cells,  the  external  emission 
source  does  not  need  to  have  a  subscript  in  order  to  define  which  cell  it  is  in. 
Both  of  the  sources  are  the  same. 

In  order  to  merge  the  two  sub-cells,  the  currents  at  the  center  face  must  be 
eliminated.  This  is  done  through  a  series  of  linear  algebra  equations.  In  this 
section  I  will  focus  on  the  merging  of  the  matrices  resulting  in  current  only.  The 
merging  of  matrices  for  the  flux  is  presented  in  Appendix  C 


(30) 

(31) 

(32) 

(33) 


Substituting  (32)  into  (31)  and  solving  for  j+c  we  get: 

1+C  =  R  +_(T .  j-^  +  R  F+j+C  +  E  .  $ext-  E  -+se^+) 

j  \  j  o  7  7  7  / 


3 


3 


_l_rp  ++~j+.Zy  H — g  ++geir£.+ 

J  j  j 


(34) 


~\*c  =  (I  -  R  II  _+r1(R  +~t  +  R  +~E  “S' 

J  v  3  3  ’  v  3  3  J  3  3 

+R  +_E  ~+sext,+ 

0,1  3 

++  •  r _j_  E  H — g ext,—  _|_  E  h — hgeait.+  N 
J  J  J  t  ' 


h — i )  — i-\— 1/t>  4 — ti  — ~l~R  ,  i >  h — i,'  — Sext,— 


(35) 


Similarly,  substituting  (31)  into  (32)  and  solving  for  j"  we  get: 
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(36) 


i“c  =  T  i~n  +  R  -+(R  +-j-u  +  T ++j+i'  +  E  +-Sm,_  +  E  ++Sm+)  + 

J  j  J  j  \  j  J  j  J  j  j  ) 


-R 


— b /T3  H — Z~c  ,  rp  ++j  +L  H — Qext,- 

— b  c\ext,+ 


E  .—Sext~  +  E  r+se 

3  3 


'fC  =  (I  -  R  ~+R  +_)_1(T.  j-^  +  R  .-+T  ++j+i  +  R  -+E  +-Sext~  + 

J  v  j  j  ’  v  j  J  j  j  J  j  j 

R  “+E  ++Sext-+  + 


3  j 
— cext,- 


- ge£t,—  _j_  g  — hg63^,+  \ 

J  J  ' 


(37) 


Substituting  (37)  into  (30)  we  get  the  current  out  the  left  side  of  the  cell: 


L  =  T  ~//j  _  R  -+R  +-)-VT“j“"  +  R  -+t  ++j+L  +  R  -+E  +~Sext~  + 

J  Jvv  3  3  K  3  J  3  3  J  3  3 


-+TD  +-\-VT - 1~R  ,  T>  -+T1  ++~\+L  ,  T>  “+T?  H  csxt,- 


R  _+E  ++sext-+  + 

3  3 

— b  Cfext,+  \\  ,  13  — b  ~i+L  ,  in QCXt,—  .  rn  — b(neiT^,+ 


E  — Sm  “  +  E  ,-+S  ’  ))  +  R  r+]+L  +  E  .  s  ,_  +  E  _+S' 

3  3  ”  3  J  3  3 


(38) 


Likewise,  substituting  (35)  into  (33)  we  get  the  current  out  of  the  right  side  of 
the  cell: 


]+R  =  r  +~yR  +  t ++((i  -  r  +  r  r+ri(  R  +-T  “j-"  +  R  +-e  r~se 

J  3  3  vv  3  3  '  v  3  3  3  3 

+R  +~E  ~+Sext+ 

J  J 

+  +  "•+£  _j_  H — g6#£,—  _j_  +  +  g6fC^.+  \\  _j_  g  H — geic£,—  _j_  ++g6X^.4 

jo  j  j  J  J  j  j 


H — T>  — b  \  1  /"T>  H — rp - ~i~R  ,  T>  H — 17 - Q\ext.— 


(39) 


Collecting  the  terms  that  operate  on  the  current  coming  in  both  sides  a  new 
m™1’"'  matrix  is  created  for  the  merged  cell  which  can  be  separated  into  new  T 
and  R  sub-matrices.  Using  the  example  from  figures  6  through  11,  using  the 
number  of  merges  n  =  3  and  going  through  the  merges  as  i  =  n  to  1  using  steps 
of  -1  the  new  sub-matrices  are: 


24 


Following  this  example,  the  original  m-’”  matrix,  at  level  n  =  3  (sub-cell 

thickness  =\  =  ~),  was 

23  8  ; 


f . 


m 


.out, in 


X(3  )j  R(3)i  + 

p  H —  rp  ++ 

V  (3 )3  X(3 )j  J 


(44) 


After  3  merges,  at  each  level  replacing  the  old  sub-matrices  with  the  new 
merged-cell  m matrix  is: 


m 


out, in 


L(0  )j 


D  — t- 

tt(o  )j 


p  H —  rp  ++ 

V  (0 )j  x(0 )j  ) 


(45) 


with  a  new  cell  thickness  of  =  1 . 

2° 

A  similar  exercise  can  be  done  for  the  terms  that  operate  on  the  external 
emission  source  going  in  each  direction  of  the  cell.  A  new  is  created  for  the 

combined  cell  that  can  also  be  separated  into  new  E  sub-matrices. 


E-  ,.r-=T„.,--(I-R„,,-"R,i,,+T1R 


(*-l  )j 


(®)i 


(i)j  “( i)j  )  “(i)j  E(*)i  + 


T/  A  ■  (I  —  R,,  A  .  +R,  A  . 

(*b  v  Ob  Ob 


+  -N-1 


E,.  “+  =  T,a  .  (I  -  R,  a  F+R  , 

(*— i)j  Ob  v  Ob  Ob 


)  E0b  +  E0b 
R/  A  ,_  +  E 


Ob 

—  i  >  H — \-ln  — Htti  ++ 


i)i  Ob 


\i)j 


(I-R.a  _+R,-)  +_)_1E/ a  ~+  +  E 

Ob  Ob  3  Ob 


Ob' 


(46) 


(47) 


E. 


i-W 


+-  -  X,  ++(I  -  Rr)  +“R,a  .“+)_1R,a  +“E, 


Ob' 


Ob 

T,a,++(I  -  R  A  +  R  •  \  •  •  )  'E  +  E  ^ 


Ob  (hi 


+ 


Ob'  R0b 


Ob 


Ob' 


(48) 


T71  ++  _  rp  ++7T  l>  H — I)  — ^ — TP  — i 

Eo-ib  ~\i)j  ^  R(i)j  Ra)j  )  R(i)j  Eob  + 


( i)j 


Ob 

T  ,++(I  -  R  ,  R  a  •  '  )  1  E  ,++  +  E  ++ 


(i)j  R0b 


Ob' 


V)j 


(49) 


Similar  to  (45)  the  merged-cell  m is: 


m 


out, ext 


E(0)J 


E<o)T 


1 71  H TT'  +  + 

v>b  E(ob  ; 


(50) 
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By  merging  the  sub-matrices  in  this  way  and  then  reassembling  the  m 
matrices  after  each  merge,  this  algorithm  can  be  repeated  n  times,  giving  a  final 
set  of  m  matrices  that  can  now  act  on  the  current  coming  into  the  cell  and  the 
external  emission  source  in  the  cell  to  produce  a  current  out  of  the  cell  in  each 
direction.  These  directions  can  then  be  summed  over  the  ordinate  weights  to 
give  a  partial  current  out  of  the  cell  in  each  direction  that  has  the  accuracy  of  211 
cells  each  of  which  is  2"n  thick.  A  simplified  version  can  be  seen  in  Algorithm  2. 
As  can  been  seen  in  the  above  matrices,  a  matrix  inversion  is  required  during 
each  merging  of  cells.  This  is  a  possible  point  of  bad  conditioning  of  the 
matrices.  By  numerically  taking  inverses  of  matrices,  loss  of  good  digits  can 
occur.  For  this  method  I  used  a  Gaussian  Elimination  with  Partial  Pivot  routine 
[16]  that  ensured  a  maximum  conditioning  of  the  matrices  during  each  merge. 
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III.  TESTING 


In  order  to  test  this  new  approach  to  DI,  it  was  necessary  to  verify  that  it 
accurately  solves  for  the  partial  currents  out  of  the  cell  and  the  average  scalar 
flux  in  the  cell.  This  is  done  by  comparing  the  solutions  from  this  algorithm  to 
solutions  from  a  benchmark  solution  achieved  from  an  SI  code.  Once  the  code  is 
shown  to  be  capable  of  achieving  the  same  solutions  as  SI  then  the  code  can  be 
tested  against  that  same  benchmark  to  demonstrate  its  ability  to  reach  the  same 
solution  with  fewer  computations  than  the  SI  code. 


III. A.  Verification  Testing 

To  verify  the  accuracy  of  my  code,  I  used  an  SI  code  [18]  that  uses  various 
input  and  material  parameters  to  solve  for  the  currents  out  of  each  cell  and  the 
average  scalar  flux  in  the  cell.  The  parameters  for  the  SI  code  include:  Ordinate 
set,  illumination  on  each  side  of  the  cell  (current  in),  emission  source  in  the  cell, 
total  cross  section  (at),  and  scattering  ratio  (c=  aj  as).  Each  of  these 
parameters  was  adjusted  in  various  combinations.  The  SI  code  was  run  until  it 
reached  a  limit  where  it  ceased  to  change  with  a  tolerance  of  10"16.  This  allowed 
me  to  compare  the  answer  between  my  DI  method  and  SI  to  16  digits.  The  SI 
code  is  capable  of  using  the  Step,  SC,  DD  and  EC  spatial  quadratures.  Because 
my  approach  is  only  concerned  with  spatial  quadratures  that  are  positive  and 
linear,  the  verification  portion  of  the  testing  only  compares  the  two  methods  that 
meet  these  criteria,  Step  and  SC.  Due  to  the  fact  that  Step  and  SC  are  only 
first-order  accurate  methods,  the  key  to  attaining  accuracy  comparable  to  the 
higher-order  methods  is  to,  in  the  case  of  SI,  create  a  finer  mesh  with  cells  that 
are  optically  thinner,  and  then  perform  the  calculations  through  the  thinner 
celled  problem.  In  section  II  of  this  document  I  presented  how  this  was  to  be 
achieved  with  my  method,  by  subdividing  the  cell,  and  then  repeatedly  merging 
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sets  of  two  sub-cells  at  a  time  until  the  original  cell  is  solved.  Various  cell 
divisions  were  tested  to  verify  that  there  was  no  loss  of  precision  in  my  method. 

Step  1  in  the  verification  process  is  to  ensure  that  the  arithmetic  of  my 
method  is  valid.  I  did  this  by  setting  up  a  series  of  conditions  in  my  method  that 
would  let  me  verify  that  the  individual  entries  in  the  m  matrices  were  correct. 

To  do  this  I  used  the  following  problem  definitions: 

•  Vacuum  boundaries  on  both  sides  of  the  cell 

•  Number  of  divisions  (for  my  method)  =  4 

•  Cell  mesh  refinement  (for  SI  code)  =  256 

•  Angular  quadrature  =  S4 

•  cr  =  .80 

•  c  =  .65 


•  At  =  20  (mfp  =  cr  At  =16) 

Using  the  matrix  form  of  (21)  and  setting  all  but  one  entry  in  jm  to  0,  and 
all  of  Sext  to  0,  j0ilt  should  equal  the  first  column  in  m0llt'in . 

'v  b  . in 

J  /ul 


Jout  _ 


(  out, in  'i 

m  1,1 

f  out, in 

m  1,1 

out, in 
m  1,2 

out, in 
m  1,3 

mouUnn 

„  out, in 
m  2,1 

out, in 

m  2,2 

mouUnx  3 

out, in 
m  3,1 

out, in 
m  3,1 

out, in 
m  3,2 

out, in 
m  3,3 

out, in 

Vil  4,1  J 

[mouUinAyX 

mout"l,\2 

mout’in43 

out, in 
m  1,4 


111 


out, in 


2,4 


out,in 
m  3,4 


ill 


out, in 


4,4  ) 


J  /u2 

■in 

J  /u3 
■in 

J  /u4 


+  m 


out, ext 


0(51) 


where  jin =  1 ,  and  jm ^ 2  =  jin ^ 3  =  jm ^ 4  =  0  .  This  process  was  then  repeated 
setting  each  consecutive  entry  of  jin  to  1  and  the  rest  to  zero,  and  then  doing  the 
same  for  the  Sext  vector.  I  used  the  same  input  variables  and  conditions  for  the 
SI  code  and  compared  the  results  (for  Step  and  SC)  for  each  column  of  the  m 
matrices.  Using  this  scheme  the  Symmetric  Relative  Difference  (SRD)  between  SI 
and  my  method  was  on  the  order  of  10"Li.  Based  on  these  results,  I  conclude  that 
my  method  is  calculating  the  within-cell  transport  correctly. 
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III.B.  Performance  Testing 

In  order  to  test  the  performance  of  the  merging  approach  to  DI,  the  EC 
method  of  SI  was  used  as  the  benchmark.  Even  though  it  is  not  a  linear  method, 
it  has  been  shown  [15]  that  it  has  a  high  order  of  accuracy  (4)  and  converges 
faster  than  any  of  the  other  methods  contained  in  the  SI  code.  It  also  has  the 
benefit  of  being  an  unconditionally  positive  quadrature. 

In  order  to  use  EC  as  a  benchmark,  I  ran  the  SI  code  multiple  times  using  an 
S8  angular  quadrature,  each  run  increasing  the  sub-cell  mesh  until  the  average 
scalar  flux  (for  test  problems  1,  2,  and  3)  and  the  partial  current  (for  test 
problem  4)  stopped  changing  at  10 (>.  For  EC  this  took  a  mesh  refinement  of  25 
(32)  sub-cells  for  a  scattering  ratio  of  c=.25  and  26  (64)  sub-cells  for  c=.5  and 
c=.75.  Once  the  benchmark  was  set,  I  repeated  the  process  for  each  of  the  other 
three  spatial  quadratures  in  the  SI  code,  ST,  SC,  and  DD,  until  the  difference 
between  it  and  the  benchmark,  for  average  scalar  flux  (test  problems  1  and  2) 
and  the  partial  current  (test  problem  3)  was  less  than  10  °.  I  repeated  this 
process  using  my  code  and  compared  the  results  of  the  SI  code  and  my  code 
(using  Step  and  SC)  in  order  to  verify  that  my  code  was  obtaining  the  same 
answers.  The  results  of  this  series  of  tests  show  that  for  several  different 
combinations  of  input  and  material  properties  that  my  method  was  consistent 
with  SI  methods  (Step  and  SC)  with  an  SRD  <  10"13  for  problems  that  had  low 
scattering  ratios  and  were  not  optically  thick. 

However,  as  the  cell  was  divided  many  times,  with  the  original  cell  thickness 
becoming  optically  thick  and  the  scattering  ratio  increasing,  as  subsequent  tests 
will  reveal,  the  answers  began  to  lose  accuracy.  Using  Step  the  two  methods 
began  to  diverge  after  7  divisions  (128  mesh  refinement)  with  an  SRD  <  10"12. 

At  10  cell  divisions  (1024  mesh  refinement)  the  two  methods  had  an  SRD  <  10"11. 
At  12  cell  divisions  (4096  mesh  refinement)  the  two  methods  had  an  SRD  <  10"10. 
Even  after  20  divisions  (1048576  mesh  refinement)  only  had  an  SRD  <  10 '.  This 
loss  of  digits  is  most  likely  due  to  calculating  the  multiple  matrix  inverses  that 
are  required  during  each  merge  of  the  m  matrices.  Even  though  the  matrices  are 
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very  well  conditioned  (matrix  condition  number  =  1  for  low  to  mid-range 
scattering  media),  the  loss  of  precision  after  7  divisions  results  in  a  negligible  loss 
of  digits.  This  is  further  investigated  in  section  III.C,  Matrix  Conditioning. 

III.B.l  Test  Problem  1:  Scattering  Ratio  (c)  =  .25 
This  test  problem  is  characteristic  of  a  high  absorbing  material  with  an 
emission  source  within  the  material.  The  parameters  for  this  test  were: 

•  Vacuum  boundaries  on  both  sides  of  the  cell 

•  Angular  quadrature  =  Sg 

•  cr  =  .75 

•  c  =  .25 

•  Ax  =  2  (mfp  =  1.5) 

•  Sext  =  2 

_  ~*-in  in  ~ 

•  Jl=Jr= 2 

The  results  of  the  test  are  contained  in  Table  3  and  a  visual  graphic  display 
shown  in  Figure  13. 

Table  2:  Test  Problem  1  Results 


Spatial  Quadrature 

Mesh  Refinement/No 
of  Merges  (No.  of  sub¬ 
cells) 

Number  of 
Iterations  (SI)/ 
Number  of 
Merges  (Dl) 

Time  to 

Reach  Limit 
[s] 

SI:  EC 

32 

23 

0.078 

SI:  DD 

64 

25 

0.063 

SI:  ST 

131072 

25 

6.45 

SI:  SC 

512 

26 

0.125 

Dl:  ST 

17 

17 

0 

Dl:  SC 

9 

9 

0.016 

As  shown  in  these  results,  for  this  type  of  material  DD  is  a  good  performer. 
However,  looking  only  at  average  scalar  flux  can  be  deceiving  in  that  the  negative 
individual  vector  components  of  the  flux  can  be  negative,  which  as  stated  before 
are  non-physical  and  hence  meaningless.  As  can  be  seen  in  the  table  Step  using 
SI  required  a  mesh  refinement  of  131,072  sub-cells  in  order  to  come  to  within  10"6 
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of  the  benchmark.  Step  using  DI,  however,  only  required  17  merges  to  achieve 
the  same  results.  This  results  in  much  less  computational  time. 


Spatial  Quadrature  <|>A  Dependence 
On  Mesh  Refinement 

S  =2,  JLinc  =JRjnc=  2,  CTt  =  .75/cm,  c  =  .25 
Ax  =  2  cm,  atAx  =  1 .5  mfp 


a) 

o 

c 

0) 

I 


cd 

> 

TO 

CD 

cr 

o 


CD 

E 

E 

>> 

C/) 


n 

SI  =  2"  sub-cells  in  mesh 
DI  =  n  sub-cell  merges 

Figure  13:  Test  Problem  1:  c  =  .25 


III.B.2  Test  Problem  2:  Scattering  Ratio  (c)  =  .75 
Using  the  same  source  and  incident  currents  as  Test  Problem  1,  this  test 
increased  the  scattering  ratio  to  0.75,  characteristic  of  a  high  scatter  medium. 
The  set  of  parameters  for  this  test  are: 

•  Vacuum  boundaries  on  both  sides  of  the  cell 

•  Angular  quadrature  =  S8 
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•  cr  =  .75 

•  c  =  .75 

•  Ax  =  2  (mfp  =  1.5) 

•  Sext  =  2 

•  jt  =  jR=2 

This  case  is  similar  to  Test  Problem  1,  but  the  average  flux  in  the  cell 
becomes  more  difficult  for  SI  to  achieve  the  correct  answer,  requiring  a  more 
refined  cell  and  more  computational  time. 


Table  3:  Test  Problem  2  Results 


Spatial  Quadrature 

Mesh  Refinement/No 
of  Merges  (No.  of  sub¬ 
cells) 

Number  of 
Iterations  (SI)/ 
Number  of 
Merges  (DI) 

Time  to  Reach 
Limit  [s] 

SI:  EC 

64 

61 

0.156 

SI:  DD 

512 

64 

0.125 

SI:  ST 

262144 

65 

20.48 

SI:  SC 

512 

67 

0.219 

DI:  ST 

18 

18 

0 

DI:  SC 

9 

9 

0 

As  shown  in  Table  4,  DD  continues  to  approach  the  same  limit  as  EC 
quickly  but  requires  twice  the  cell  refinement  in  order  to  do  so.  SC  in  this  case 
requires  no  more  refinement  to  approach  the  same  limit  as  EC  as  before  in  test  1, 
but  does  require  more  computational  time  because  of  the  increase  scattering 
ratio.  In  both  cases  of  the  DI  routines,  even  though  the  number  of  merges 
increased,  the  effective  computational  cost  did  not.  This  shows  the  strength  of 
this  new  approach.  Figure  14  shows  a  similar  pattern  for  this  test  as  for  test  1. 
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Spatial  Quadrature  <jrx  Dependence 
On  Mesh  Refinement 
S  =2,  JLinc  =JRinc=  2,  Ct  =  .75/cm,  c  =  .75 
Ax  =  2  cm,  atAx  =  1 .5  mfp 


Figure  14:  Test  Problem  2:  c  =  .75 

III.B.3  Test  Problem  3:  No  Source,  Current  in  on  One  Side 
In  this  test,  the  parameters  of  the  problem  were  changed  to  be  a  medium  with 
no  external  source,  incident  current  on  the  left  side  of  the  cell  only,  and  the 
width  of  the  cell  beginning  at  3  mfp  instead  of  1.5.  The  results  of  this  test  are 
summarized  in  Table  5.  The  set  of  parameters  for  this  test  are: 

•  Vacuum  boundaries  on  both  sides  of  the  cell 

•  Angular  quadrature  =  S8 

•  cr  =  .75 

•  c  =  .75 
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•  Ax  =  4  (mfp  =  3) 

•  Sext  =  0 

•  Tl=  2 

•  7f  =  o 


Table  4:  Test  Problem  3  Results 


Spatial  Quadrature 

Mesh  Refinement/No 
of  Merges  (No.  of  sub¬ 
cells) 

Number  of  Iterations 
(Sl)/Number  of 
Merges  (Dl) 

Time  to  Reach  Limit 
[s] 

SI:  EC 

64 

88 

0.359 

SI:  DD 

1024 

93 

0.594 

SI:  ST 

262144 

89 

33.4 

SI:  SC 

512 

89 

0.594 

Dl:  ST 

18 

18 

0 

Dl:  SC 

9 

9 

0 

The  results  of  this  test  show  the  weakness  of  DD  in  that  when  the  cell  is 
optically  thick  (>2  mfp)  the  partial  currents  on  the  far  side  of  the  cell  are 
negative  until  the  mesh  refinement  is  high  enough.  This  negative  current  is 
shown  in  Figure  15. 
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Spatial  Quadrature  jOUI  Dependence 
On  Mesh  Refinement 

S  =0,  JLinc  =  2,  JRinc=  0,  a,  =  .75/cm,  c  =  .75 
Ax  =  4  cm,  o', Ax  =  3  mfp 


0  2  4  6  8  10  12 

n 

SI  =  2n  sub-cells  in  mesh 
Dl  =  n  sub-cell  merges 

Figure  15:  Test  Problem  3:  No  Soure,  One  Current 


III.B.f  Test  Problem,  4:  Approaching  Diffusion 

The  final  test  that  was  conducted  is  one  that  tests  the  limits  of  discrete 
ordinates  methods,  that  of  a  diffusion  problem  where  c  =  f  and  the  cell  is  many 
mean  free  paths  thick.  In  order  to  do  this,  I  set  the  problem  parameters  to: 

•  Vacuum  boundaries  on  both  sides  of  the  cell 

•  Angular  quadrature  =  S8 

•  cr  =  .75 

•  c  =  .99 

•  Ax  =  2  (mfp  =  1.5) 

•  Sext  =  0 
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_  ~*'in  “? in  r\ 

•  Jl=Jr=  0 

The  benchmark  SI  code  used  was  not  able  to  reach  a  limit  for  c  =  1. 
Therefore,  for  this  test,  the  scattering  ratio  only  approaches  1,  hence  it 
approaches  diffusion.  Even  at  c  =  .99  only  the  SI  quadrature  of  EC  was  able  to 
reach  a  limit  in  a  timely  manner  requiring  over  1100  sweeps  to  stop  changing  at 
10"G.  This  brings  to  question  the  validity  of  the  results,  but  serves  as  a  decent 
comparison  between  the  computational  speeds  of  this  method  and  those  of  SI. 
The  results  of  this  test  are  shown  in  Figure  16. 

Partial  Current  Convergence 
Diffusion  Problem 

S  =  2,  JLjnc  =  0,  JRjnc  =  0,  at  =  .75/cm,  c  =  .99 
Ax  =  100  cm,  ctjAx  =  75  mfp 


SI  =  2n  sub-cells  in  mesh 
Dl  =  n  sub-cell  merges 


Figure  16:  Test  Problem  4:  Approaching  Diffusion 


36 


As  can  be  seen  from  the  figure,  DI:SC  reached  the  same  limit  for  partial 
current  as  SI:EC.  As  stated  above,  EC  required  over  1100  sweeps  to  reach  this 
limit,  and  this  had  to  be  done  through  210  sub-cells.  This  required  a  real  time 
solution  of  approximately  1  minute.  DI  required  13  merges  to  reach  the  same 
solution,  requiring  computational  time  of  0.0156  seconds. 

III.C.  Matrix  Conditioning 

As  stated  earlier,  the  number  of  matrix  inversions  is  a  possible  point  of  loss  of 
precision  in  my  method.  The  following  tests  show  how  the  condition  number  of 
the  matrix  varies  with  different  problem  parameters,  including  the  number  of 
merges  required  to  achieve  the  accuracy  of  high-order  spatial  quadratures.  In 
each  of  these  tests  I  used  the  infinity-norms  of  the  matrices  in  order  to  calculate 
the  condition  numbers. 

III.C.  1  Condition  Number  vs.  Cell  Thickness  (Optical)  (c  =  .1) 

The  first  parameter  to  investigate  is  cell  thickness  in  terms  of  the  optical 
thickness  of  the  medium.  The  parameters  for  this  set  of  tests  are  the  following: 

•  Vacuum  boundaries  on  both  sides  of  the  cell 

•  Angular  quadrature  =  S8 

•  cr  =  1 

•  c  =  .1 

•  Original  cell  At  (varies)  =  1  to  1024  mfp 

•  Sext  =  5 

•  Jr = 2 

•  7*  =3 

•  Number  of  merges  =  20 

Based  on  the  results  of  this  test  the  condition  number  went  no  higher  than 
1.0017.  This  demonstrates  that  a  problem  with  a  low  scattering  ratio  the  cells 
can  be  very  optically  thick  without  increasing  the  condition  number  with  20 
merges.  My  method  is  well  conditioned  for  this  set  of  parameters.  For  this 
problem,  using  a  high  absorbing  medium,  my  method  can  take  a  cell  that  has  an 
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optical  thickness  of  1024  mean  free  paths,  and  using  a  first-order,  linear,  positive 
spatial  quadrature,  achieve  the  same  accuracy  as  if  the  cell  was  1/1000  the 
thickness  of  a  mean  free  path. 

III.C.2  Condition  Number  vs.  Cell  Thickness  (Optical)  (c  =  .9) 

In  this  test  I  use  the  same  parameters  as  those  in  the  last  section  except  for 
the  scattering  ratio,  which  I  now  set  at  c  =  .9.  This  is  a  high  scattering  medium. 
The  results  for  this  test  are  shown  in  Table  6.  In  this  table  the  first  column  is 
the  merge  number.  For  example,  merge  #1  is  the  first  merge  of  two  of  the  most 
optically  thin  sub-cells.  Merge  ((‘20  is  the  last  merge  of  the  two  half  sub-cells  to 
make  result  in  the  original  cell.  This  is  true  for  both  of  the  next  two  tables. 

The  results  for  this  test  are  comparable  to  the  previous  test,  showing  that 
within  the  range  of  scattering  ratios,  from  c  =  .1  to  c  =  .9,  the  thickness  of  the 
cell  does  not  have  a  significant  effect  on  the  conditioning  of  the  matrices. 
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Table  5:  Condition  Number  vs.  Cell  Thickness  (c  =  0.9) 


Cell  Thickness  (Optical)  c  =  .9 

Merge 

# 

1 

2 

4 

8 

16 

32 

64 

128 

256 

512 

1024 

1 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

2 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

3 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0001 

4 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0001 

1.0002 

5 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0001 

1.0002 

1.0008 

6 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0001 

1.0002 

1.0008 

1.0029 

7 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0001 

1.0002 

1.0008 

1.0029 

1.0105 

8 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0001 

1.0002 

1.0008 

1.0029 

1.0105 

1.0342 

9 

1.0000 

1.0000 

1.0000 

1.0000 

1.0001 

1.0002 

1.0008 

1.0029 

1.0105 

1.0342 

1.0980 

10 

1.0000 

1.0000 

1.0000 

1.0001 

1.0002 

1.0008 

1.0029 

1.0105 

1.0343 

1.0982 

1.2417 

11 

1.0000 

1.0000 

1.0001 

1.0002 

1.0008 

1.0029 

1.0105 

1.0343 

1.0982 

1.2419 

1.5086 

12 

1.0000 

1.0001 

1.0002 

1.0008 

1.0029 

1.0105 

1.0343 

1.0982 

1.2420 

1.5089 

1.8572 

13 

1.0001 

1.0002 

1.0008 

1.0029 

1.0105 

1.0343 

1.0983 

1.2421 

1.5091 

1.8577 

2.0659 

14 

1.0002 

1.0008 

1.0029 

1.0105 

1.0343 

1.0983 

1.2421 

1.5091 

1.8579 

2.0663 

2.0974 

15 

1.0008 

1.0029 

1.0105 

1.0343 

1.0983 

1.2421 

1.5092 

1.8580 

2.0665 

2.0978 

2.0979 

16 

1.0029 

1.0105 

1.0343 

1.0983 

1.2421 

1.5092 

1.8580 

2.0667 

2.0980 

2.0983 

2.0979 

17 

1.0105 

1.0343 

1.0983 

1.2421 

1.5092 

1.8581 

2.0667 

2.0981 

2.0985 

2.0983 

2.0979 

18 

1.0343 

1.0983 

1.2421 

1.5092 

1.8581 

2.0667 

2.0982 

2.0986 

2.0985 

2.0983 

2.0979 

19 

1.0983 

1.2421 

1.5092 

1.8581 

2.0668 

2.0982 

2.0987 

2.0986 

2.0985 

2.0983 

2.0979 

20 

1.2421 

1.5092 

1.8581 

2.0668 

2.0982 

2.0987 

2.0987 

2.0986 

2.0985 

2.0983 

2.0979 
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III.C.3  Condition  Number  vs.  Cell  Thickness  (Optical)  (c  =  1) 

This  test  focuses  on  a  pure  scattering  medium,  a  problem  that  is  not  possible 
for  SI  codes  to  solve.  Using  all  the  same  parameters,  changing  the  scattering 
ratio  to  c  =  1,  I  ran  the  same  problem,  using  cell  thicknesses  from  1  mean  free 
path  to  1024  mean  free  paths.  The  results  are  shown  in  Table  7. 


Table  6:  Condition  Number  vs.  Cell  Thickness  (c  =  1) 


Cell  Thickness  (Optical)  c  =  1 

Merge 

# 

1 

2 

4 

8 

16 

32 

64 

128 

256 

512 

1024 

1 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

2 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

3 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0001 

4 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0001 

1.0002 

5 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0001 

1.0002 

1.0010 

6 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0001 

1.0002 

1.0010 

1.0037 

7 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0001 

1.0002 

1.0010 

1.0037 

1.0132 

8 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0001 

1.0002 

1.0010 

1.0037 

1.0132 

1.0440 

9 

1.0000 

1.0000 

1.0000 

1.0000 

1.0001 

1.0002 

1.0010 

1.0037 

1.0132 

1.0440 

1.1320 

10 

1.0000 

1.0000 

1.0000 

1.0001 

1.0002 

1.0010 

1.0037 

1.0132 

1.0441 

1.1321 

1.3570 

11 

1.0000 

1.0000 

1.0001 

1.0002 

1.0010 

1.0037 

1.0132 

1.0441 

1.1322 

1.3573 

1.9107 

12 

1.0000 

1.0001 

1.0003 

1.0010 

1.0037 

1.0132 

1.0441 

1.1322 

1.3574 

1.9113 

3.2696 

13 

1.0001 

1.0003 

1.0010 

1.0037 

1.0132 

1.0441 

1.1322 

1.3575 

1.9116 

3.2707 

6.4441 

14 

1.0003 

1.0010 

1.0037 

1.0132 

1.0441 

1.1322 

1.3575 

1.9117 

3.2713 

6.4466 

13.3319 

15 

1.0010 

1.0037 

1.0132 

1.0441 

1.1322 

1.3575 

1.9118 

3.2716 

6.4479 

13.3372 

27.5573 

16 

1.0037 

1.0132 

1.0441 

1.1322 

1.3575 

1.9118 

3.2718 

6.4485 

13.3398 

27.5678 

56.3169 

17 

1.0132 

1.0441 

1.1322 

1.3576 

1.9118 

3.2718 

6.4488 

13.3411 

27.5731 

56.3380 

114.0213 

18 

1.0441 

1.1322 

1.3576 

1.9118 

3.2719 

6.4490 

13.3418 

27.5758 

56.3486 

114.0636 

229.5322 

19 

1.1322 

1.3576 

1.9118 

3.2719 

6.4491 

13.3421 

27.5771 

56.3539 

114.0848 

229.6167 

460.6076 

20 

1.3576 

1.9118 

3.2719 

6.4491 

13.3423 

27.5778 

56.3566 

114.0954 

229.6590 

460.7765 

922.7858 
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The  results  of  this  test  show  that  for  a  scattering  medium  the  matrices 
become  ill-conditioned  as  the  cell  thickness  increases  and  the  number  of  merges 
increases.  I  shaded  a  rough  diagonal  in  Table  7,  showing  a  possible  acceptable 
limit  of  my  method  with  these  parameters  where  the  matrices  are  well- 
conditioned.  This  shows  that  for  a  given  cell  thickness  a  certain  number  of 
merges  can  be  performed  before  an  unacceptable  loss  of  good  digits  occurs.  For 
example,  if  I  were  to  begin  with  a  cell  of  optical  thickness  =  256  mean  free  paths, 
I  could  expect  to  be  able  to  merge  sub-cells  that  are  2"14  -  2"1'1  the  thickness  of  the 
original  cell  before  the  results  are  suspect. 

III.C.4  Condition  Number  Dependence  on  Angular  Quadrature 
The  dependence  of  the  matrix  condition  number  on  the  angular  quadrature 
used  is  a  straightforward  exercise.  As  the  number  of  ordinates  increases,  the  size 
of  the  matrix  increases.  For  example,  S4  creates  a  4  x  4  m  matrix.  Sg  creates  a  8 
x  8  m  matrix.  Because  the  m  matrices  are  almost  diagonal,  with  values  near  1 
011  the  diagonal  and  very  small  numbers  off  the  diagonal,  the  condition  number 
does  not  increase  and  is  therefore  not  dependant  on  the  number  of  ordinates,  for 
low  to  mid-range  scattering  media.  However,  as  seen  in  Table  7,  when  I  tested  a 
pure  scattering  medium  the  condition  number  increased  with  the  number  of 
merges.  After  testing  S4  -  S32  I  found  that  the  condition  number  is  not 
dependant  on  the  number  of  ordinates,  but  rather  on  the  scattering  ratio  of  the 
medium.  The  following  is  a  list  of  approximate  condition  numbers  for  the  low- 
scattering  media  through  all  20  merges: 

Angular  Quadrature  Approximate  Condition  Number 

S4  1 

S8  1 

S16  1 

S32  1 
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IV.  SUMMARY 


IV. A.  Achievement  of  Objectives 

Increased  Accuracy:  As  has  been  shown,  in  order  to  increase  accuracy  for 
the  linear,  positive  spatial  quadratures  the  thickness  of  the  cell  has  to  be 
decreased.  In  the  past,  for  SI  models,  this  has  been  done  by  refining  the  spatial 
mesh  into  many  smaller  cells.  This  increases  the  computational  time  required  to 
sweep  through  the  mesh.  For  Slab  geometry  the  increase  can  be  very  substantial. 
Extending  that  to  2  and  3-D  geometries  can  become  devastating. 

This  new  algorithm  for  DI  shows  that  it  can  achieve  the  same  accuracy 
using  Step  and  SC  spatial  quadratures  as  that  of  a  higher  order  method.  The 
bonus  to  this  fact  is  that  it  required  orders  of  magnitude  fewer  computations  to 
achieve  the  same  result. 

Numerical  Stability:  Going  into  this  project,  knowing  that  the  algorithm 
required  matrix  inversions  for  each  merge  process,  I  knew  that  stability  issues 
might  arise.  Matrix  inversions,  if  not  done  carefully,  can  result  in  the  loss  of 
digits.  For  the  range  of  optical  thicknesses  and  angular  quadrature  sets  with 
isotropic  scattering  tested  here,  the  matrices  were  remarkably  well-conditioned 
resulting  in  very  little  loss  of  good  digits  in  the  number  of  merges  performed  to 
attain  the  required  accuracy,  except  where  the  medium  was  very  optically  thick 
and  diffusive. 

Integration  into  DI  Code:  The  method  that  I  developed  has  not  been 
integrated  with  the  existing  1-D  DI  code.  The  DI  code  computes  the  mou  ,in 
matrix  before  solving  the  global  partial  current  problem  which  my  method  could 
be  used  in  place  of  the  current  method.  It  does  not,  however,  treat  the 
coefficients  of  the  external  emission  source  as  matrix  in  the  same  way  that  I  did. 
Generating  both  matrices  was  a  key  part  of  the  merging  process.  This  may 
require  some  modification  of  the  DI  code  in  order  for  it  to  use  my  method. 
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IV.  B.  Future  Work 


As  discussed  this  method  requires  integration  into  the  1-D  DI  code. 

Extending  this  method  to  2  and  3-D  geometries  will  require  much  more  work,  but 
if  achievable,  can  produce  a  new  method  that  will  out-perform  SI  codes  that  are 
currently  being  used.  Conditioning  of  the  matrices  can  be  improved  by  applying 
different  matrix  inversion  schemes.  This  will  allow  more  cell  merges,  if  necessary, 
allowing  the  starting  cell  thickness  to  be  more  mean  free  paths  thick. 


IV.  C.  Observations  and  Conclusions 
In  conclusion,  this  new  algorithm  that  can  be  used  by  the  DI  discrete 
ordinates  method  does  show  promise  in  its  ability  to  increase  the  accuracy  of 
spatial  quadratures  that  are  linear  and  unconditionally  positive.  This  will 
overcome  the  some  of  the  issues  that  discrete  ordinates  methods  have  had, 
namely  non-physical  negative  results  that  come  from  spatial  quadratures  that  are 
not  unconditionally  positive  and  solutions  that  do  not  approach  an  appropriate 
diffusion  limit  that  happens  with  non-linear  methods.  The  results  showed  very 
good  numerical  stability  and  very  short  computational  times. 
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APPENDIX  A 


K  MATRICES  FOR  STEP 

The  K  matrices  are  defined  by  the  spatial  quadrature  and  the  balance  equation. 
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cr  Ax 
A, i 
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As  is  explained  by  Mathews  [14],  for  an  isotropic  source,  like  the  source  in  the 
examples  that  I  used  to  test  my  method  the  following  applies: 
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All  of  these  matrices  are  diagonal  and  symmetric. 
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APPENDIX  B 

K  MATRICES  FOR  STEP  CHARACTERISTIC 
The  K  matrices  are  defined  by  the  spatial  quadrature  and  the  balance  equation. 
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where  the  zeroth  exponential  moment  function  is  defined  as: 
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where  the  first  exponential  moment  function  is  defined  as: 
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Because  the  problem  is  isotropic  SC  follows  Step  in  that: 
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APPENDIX  C 
M  MATRICES  FOR  FLUX 

The  algorithms  for  computing  and  merging  the  m  coefficient  matrices  that  act 


on  flux  are  similar  to  those  that  were  presented  in  Section  II.C.l  of  the  thesis. 
However,  because  flux  coming  into  the  cell  creates  flux  and  current,  the  matrices 
for  flux  are  slightly  more  complicated. 

Beginning  with  equations  (25)  and  (26)  I  split  the  cell  into  two  identical  sub¬ 
cells  of  V2  dx  thicknesses  and  have  these  two  equations  for  each  cell  with  currents 
at  the  center: 
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The  flux  that  will  be  created  by  merging  the  cells  will  be  the  average  flux  across 
the  new  cell  (L  +  R)/2.  To  begin  the  process  I  add  (47)  to  (49)  and  (48)  to  (50) 
which  gives  the  following: 
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and: 
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Substituting  (35)  and  (37)  into  both  of  these  equations  in  order  to  eliminate  the 
currents  at  the  center  boundary  gives  the  following: 
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I  repeat  the  procedure  outlined  in  Section  II.C.l  for  the  current  equations, 
collect  the  terms  that  act  on  the  current  coming  in  both  sides  of  the  cell,  as  well 
as  the  flux  and  rearrange  to  get  new  m  matrices.  For  the  mv/,in  matrix  that  acts 
on  the  current  I  get: 
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Now  the  new  mv,in  matrix  is: 
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Doing  the  same  for  the  external  emission  source  terms  yields: 
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Now  the  new  merged-cell  m  matrix  for  flux  due  to  external  emission  source  is: 
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